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ABSTRACT 

This paper uses dynamical invariants to describe the evolution of collisionless systems subject 
to time-dependent gravitational forces without resorting to maximum-entropy probabilities. 
We show that collisionless relaxation can be viewed as a special type of diffusion process 
in the integral-of-motion space. In time-varying potentials with a fixed spatial symmetry the 
diffusion coefficients are closely related to virial quantities, such as the specific moment of 
inertia, the virial factor and the mean kinetic and potential energy of microcanonical particle 
ensembles. The non-equilibrium distribution function (DF) is found by convolving the ini¬ 
tial DF with the Green function that solves Einstein’s equation for freely diffusing particles. 
Such a convolution also yields a natural solution to the Fokker-Planck equations in the en¬ 
ergy space. Our mathematical formalism can be generalized to potentials with a time-varying 
symmetry, where diffusion extends over multiple dimensions of the integral-of-motion space. 
The new probability theory is in many ways analogous to stochastic calculus, with two signif¬ 
icant differences: (i) the equations of motion that govern the trajectories of particles are fully 
deterministic, and (ii) the diffusion coefficients can be derived self-consistently from micro- 
canonical phase-space averages without relying on ergodicity assumptions. For illustration we 
follow the cold collapse of A^-body models in a time-dependent logarithmic potential. Com¬ 
parison between the analytical and numerical results shows excellent agreement in regions 
where the potential evolution does not depart too strongly from the adiabatic regime. 

Key words: Galaxy: kinematics and dynamics; galaxies: evolution; diffusion < Physical Data 
and Processes. 


1 INTRODUCTION 

Systems with a very large number of degrees of freedom exhibit 
regularities that can be expressed as statistical laws. For gases and 
plasmas it is possible to derive such laws in a reasonably straight¬ 
forward manner from certain general principles of statistical me¬ 
chanics (e.g. Landau & Lifshitz 1980). Unfortunately, statistical 
mechanics has serious, unresolved issues when it comes to de¬ 
scribing the non-equilibrium evolution of dynamical systems sub¬ 
ject to long-range forces such as gravity (see Padmanabhan 1990 
and Campa et al. 2009 for detailed reviews). Particles interacting 
via long-range forces exhibit some unusual behaviour, like non- 
ergodicity and negative specific heat (Antonov 1961; Lynden-Bell 
& Lynden-Bell 1977; see Lynden-Bell 1999 for a review), which 
is difficult to incorporate into a classical thermodynamical frame¬ 
work. These difficulties are in part responsible for the current lack 
of theoretical understanding of how virialized systems settle in their 
particular equilibrium configuration. 

The first attempt to construct a statistical theory for the evo¬ 
lution of collisionless galaxies is due to Lynden-Bell (1967), who 
used concepts of classical thermodynamics in order to predict the^- 
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nal equilibrium state of a collisionless gravitating system evolving 
in a time-dependent gravitational field, a process dubbed as ‘vio¬ 
lent relaxation’. In his pioneer work Lynden-Bell discretizes phase 
space using elements (macrocells) of equal volume and different 
masses. Applying standard techniques of mechanical statistics he 
shows that violent relaxation leads toward a unique coarse-grained 
distribution function (fc) which can be described by a weighted 
superposition of Fermi-Dirac distributions with different tempera¬ 
tures. Nakamura (2000) proposes an alternative approach where he 
derives the equilibrium distribution function of macrocells of equal 
mass and different phase-space volumes using Jaynes (1957) in¬ 
formation theory. The new phase-space partition leads to a coarse¬ 
grained distribution function that follows the well-known Mawell- 
Boltzmann distribution. To date it remains unclear which of the two 
approaches is more correct from a statistical point of view (Arad & 
Lynden-Bell 2005), although both fail to reproduce numerical ex¬ 
periments of violent relaxation (Arad & Johansson 2005). 

The central difficulty in accepting the distribution function de¬ 
rived from Lynden-Bell and Nakamura’s theories is that it predicts 
infinite mass for the system. In other words, the variational prob¬ 
lem that determines the most probable fc possesses no solution for 
any finite total mass. This shortcoming may be due to the short life 
of the process that drives relaxation, i.e. fluctuations of the gravi- 
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tational field, which vanish on the time scale well before 

the thermodynamical equilibrium is attained. As a result, in most 
gravitating systems the evolution will be frozen in a subdomain of 
the available phase space (a.k.a. ‘incomplete relaxation’), break¬ 
ing the ergodicity hypothesis, which establishes an equivalence be¬ 
tween time average and average over particle ensembles. As a pos¬ 
sible remedy, Robert & Sommeria (1992) (see also Chavanis et al. 
1996) propose the ‘maximum entropy production principle’, which 
assumes that the evolution a system proceeds in such a way that 
the rate of entropy production is maximized while satisfying all the 
dynamical constraints. In this theory the coarse-grained distribu¬ 
tion function is a solution to a modified version of the (collision¬ 
less) Fokker-Planck equation. These authors find that fc evolves to¬ 
ward a Fermi-Dirac distribution as the system approaches complete 
relaxation. More recently, there have been attempts to derive the 
most probable, maximum-entropy fc after imposing additional con¬ 
straints on the system, e.g. adiabatic invariance of actions (Pontzen 
& Governato 2013). 

Unfortunately, the concept of thermodynamical ‘entropy’ {Tf) 
is ill defined in systems subject to unshielded, long range forces. 
For ideal gases the Tf-function must obey very special proper¬ 
ties (see Appendix A of Jaynes 1957) which can only be met 
by the Boltzmann entropy, T-L = -/ /cln/cd^rd^v. In contrast, for 
collisionless gravitating systems Tremaine et al. (1986) demon¬ 
strate that any convex function C'(fc) will lead to an entropy H = 
-J CM^rd^v that increases with time. In this regard the Boltzmann 
entropy C' = /cln/c is only one of a wide variety of -functions 
that must increase monotonically during violent relaxation. 

This conclusion has been rejected by Dejonghe (1987) and 
Sridhar (1987) who provide counter-examples of systems in which 
77 oscillates during mixing. Dejonghe (1987) concludes that 
Tremaine et al. theorem only proves that the 77-functions increase 
as a result of coarse-graining (which has an effect akin to infor¬ 
mation loss), and that whether or not the variation of 77 happens 
monotonically depends on (i) the prior information on the state 
of the system (e.g. the assignment of a priori probabilities to the 
macrocells) and (ii) the functional form of C'(fc). Since no partic¬ 
ular choice can be justified from variational principles it remains 
unclear how to interpret thermodynamical entropies. 

Padmanabhan (1990) identifies a different, but very worri¬ 
some limitation of thermodynamical approaches to violent relax¬ 
ation. While in systems interacting via short-range forces the en¬ 
ergy will be an extensive parameter (i.e. the total energy of the sys¬ 
tem will be the sum of the energies of the subsystems), in gravitat¬ 
ing objects this cannot be guaranteed. Indeed, the extensive nature 
of the energy does not hold in the presence of long range interac¬ 
tions because the energy of a particle in a given macrocell cannot 
be shielded from the gravitational attraction of particles in neigh¬ 
bour cells. When this happens the system can no longer be divided 
into non-interacting subsystems and the laws of equilibrium ther¬ 
modynamics do not apply. Recent work by Levin et al. (2008, 2014) 
supports this conclusion. With the aid of A-body simulations these 
authors follow the evolution of self-gravitating models with ‘water- 
bag’ distribution functions that are not initially in balance. The N- 
body models are left to oscillate with decreasing amplitude until 
they reach dynamical equilibrium. Levin et al. (see also Joyce et al. 
2009 and Sylos Labini et al. 2015) show that the final state after 
violent relaxation is typically characterized by a bi-modal distribu¬ 
tion function, with particles in the inner-most region of the poten¬ 
tial distributing according to Lynden-BelTs theory while those in 
the outer-most regions forming a distinct dilute halo which forms 


through resonances arising during the macroscopic oscillations of 
the system. 

This contribution seeks to build a probability theory for non¬ 
equilibrium gravitating systems which does not rely on thermo¬ 
dynamical arguments (e.g. phase-space discretization, maximum- 
entropy probabilities,...) or on ergodicity assumptions. Our goal is 
not only to describe the final state of violent relaxation but also 
all the intermediate stages (i.e. incomplete relaxation). We shall 
work in the collisionless, mean field limit, where the granularity 
in the system may be ignored, such that the (fine-grained) distri¬ 
bution function /(r,v,f) is assumed to be smooth over the mean 
interparticle distance of the system. 

The method proposed below follows up the work of Penarru¬ 
bia (2013, hereinafter P13), who devises of a general coordinate 
transformation that removes the explicit time-dependence from the 
equations of motion of a collisionless system (see §2). By refer¬ 
ring the phase-space locations of particles to a coordinate frame 
in which the potential remains ‘static’ the dynamical effects intro¬ 
duced by the time evolution vanish. In §3 we apply similar argu¬ 
ments to derive an invariant distribution function from the initial 
state of the system at f = to- Transforming the invariant function 
back to the original integral-of-motion space yields the desired, 
non-equilibrium /(r,v,to-l-r'), where r' > 0 is a time interval of 
arbitrary length. Interestingly, we find that this transformation is 
analogous to a Green convolution of the invariant function with 
the solution to a diffusion equation, suggesting that collisionless 
relaxation can be viewed as a special type of diffusion process in 
the integral-of-motion space. In §4 we run a suite of A-body ex¬ 
periments which follows the evolution of tracer particles under¬ 
going cold collapse in a time-dependent scale-free potential. We 
will show that the theory works well in regions where the potential 
evolution does not depart too strongly from the adiabatic regime. 
The great virtue of scale-free potentials is that the diffusion co¬ 
efficients can be expressed analytically in terms of virial quantities 
associated to microcanonical ensembles. Indeed, it will become ev¬ 
ident that deriving the diffusion coefficients in arbitrary potentials 
is generally a very difficult task. In §5 we discuss future prospects, 
like the analysis of self-gravitating systems with a time-varying 
symmetry and the possibility to incorporate the effects of random 
particle-particle encounters into our theory. For clarity Appendix A 
describes a number of concepts encountered in statistical mechan¬ 
ics, such as microcanonical ensembles, phase mixing, ergodicity,.., 
which appear repeatedly in the present work. Appendix B shows 
that Markov chains provide useful statistical tools to describe the 
evolution of collisionless systems. 


2 FROM MICRO TO MACRO 

Consider a microcanonical ensemble of particles that at the time 
t = to move on orbits with an energy £" = Fo in a time-dependent 
potential ‘I>(r,f). In this Section we use dynamical invariants (i.e. 
quantities that are conserved along the phase-space path of a parti¬ 
cle) to demonstrate that the evolution of the microcanonical distri¬ 
bution can be described by the same phenomenological equations 
that arise in stochastic processes such as Brownian motion and ran¬ 
dom walks (see Wax 1958 for a review). 
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2.1 Dynamical invariants 

The equations of motion of a particle subject to a time-dependent 
gravitational acceleration can be written as 

r = F(r,f), (1) 


where F is a specific force. If the force is conservative one can find 
a gravitational potential <1? such that F(r,f) = 

Lynden-Bell (1982) and P13 show that under these conditions 
one can construct a canonical transformation r r' R{t), and a 
time-coordinate transformation dt i—> drf?^(f), such that the explicit 
time-dependence of the force field vanishes from the equations of 
motion. In the new coordinates Equation (1) becomes 

g=r'[r'l, (2, 

Equating (1) and (2) shows that the scaling factor R{t) must be 
a solution of the following differential equation 

f?'/?V'-f?'*F(f;r',0 = -F'[r']- (3) 


Note that if F is a conservative force (V x F = 0) then F' is also con¬ 
servative. This calls for the definition of a time-independent scalar 
potential, = -f F'dr', such that in the transformed coordinates 
the energy, 1=1 /2(dr'/dr)^ + $^(r'), becomes an exact dynamical 
invariant (i.e. a constant of motion). Expressing the energy invari¬ 
ant 1 in terms of the original coordinates yields 

I = ]^{Ry-Rrf+]^RR/+R^^{r,ty, (4) 

where v = dr/df. It is important to emphasize that the energy / is 
conserved along the phase-space path of a particle motion, i.e. the 
subset of the 6-dimensional phase space on which the equations of 
motion (1) are fulfilled. 

Equation (4) can be easily re-arranged to express the energy 
as a function of the invariant quantity, that is 


E = Ea + {r-y) 


R 

R 




(5) 


where Ea(t) = I/R^(t) denotes the adiabatic energy. In what fol¬ 
lows it is useful to define the quantity A(r,v,0 = E-E^ as the 
change of energy with respect to the adiabatic solution. In gravita¬ 
tional fields that do not vary rapidly with time the energy change 
is |A/£a| < 1, which defines the so-called adiabatic regime. Note 
that A can be positive or negative depending on the sign of r • v, 
and that the tangential velocity component does not contribute to 
A, which implies that the energy change is independent of the an¬ 
gular momentum of the orbit. 


2.1.1 Scale-free potentials 

Power-law gravitational fields are of special interest in this work, as 
they permit simple approximate solutions to Equation (3). Consider 
the time-dependent power-law force 

E{r,t) = -p{t)r". (6) 


The potential associated with Equation (6) can be written as 


<E’(Tf)-‘Foo 


J n+l 

\p.(t)ln(r) 


,n^-l 
,n = -l, 


(7) 


Through a straightforward manipulation of Equation (3) P13 
shows that the scale factor 


R{t) = 


pif) 

Po 


-l/(n+3) 


(8) 


provides a solution of Equation (3) that is accurate at order Cl(e'^), 
where po = p(lo), ^' = if/ UolP and P is the radial period of the orbit 
at f = to. From Equation (8) it follows that f?//? = -(p/p)/(n-\-3). 


2.2 The microcanonical ensemble 


Let us now consider an ensemble of particles with energy E = Ea 
at f = fo. A derivation of the energy distribution at an arbitrary time 
t = to-i-r' may appear as an impossibly difficult task given that the 
value of A = £■-£„ varies according to the phase-space location of 
each individual particle in the ensemble. To attack this problem we 
shall follow a statistical method originally introduced by Einstein 
(1905)*, which is commonly applied to a wide range of stochas¬ 
tic processes encountered in Physics, Biology and Chemistry. This 
approach requires that the evolution of the gravitational field hap¬ 
pens in the linear regime and does not deviate too strongly from the 
adiabatic limit. 

As a first step one defines (/sfAliiaidA as the (conditional) 
probability that during a time interval (to,fo + 'r*) a particle with en¬ 
ergy E„ will experience an energy change A within the range dA. 
The function p is known as the probability of energy change and it 
is normalized such that 

j p{A\Ea)A/S=l. (9) 

Since A = A(r, v, t) the probability function p must be intimately 
related to the distribution of particles in phase space, an issue that 
is discussed in detail in Section 2.3. 

As a second step one introduces the function p(E,t\Ea,to), 
which defines the probability that a particle with energy Eaditt = to 
will have an energy in the range E,E-\-dE at the time t = to-tr'. 
Since the number of particles is conserved these two functions must 
obey the following equality 

p(E,to-tT'\Ea,to)dE = dE f p(E-A,to\Ea,to)p{A\Ea)dA, (10) 


which is typically known as Einstein’s master equation. 

If the evolution of the potential does not deviate too strongly 
from the adiabatic limit, i.e. | A/iia| < 1, one can expand the right- 
hand side of Equation (10) in series of A, which after taking into 
account the normalization (9) yields 


p(E,to + T'\Ea,to) 


p(E,to\Ea,to)- 


dp 


dE 

1 d^p 

2d^ 


I (p(A|£„)AdA (11) 
p(A\Ea)A^dA. 


Given that the time interval can be arbitrary small, the left-hand 
side of Equation (10) can be also approximated as 

p{E,to+T'\Ea,to) ~ p(E,to\Ea,to)-l-^T' , (12) 

Equating (11) and (12) yields 


dp^ d[^ d^p 

dt dE dE^ 


(13) 


An English translation can be found in Stachel (1998). 


where ^oo = const. 


1 
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where the drift coefficient is defined as 

C{Ea,t) = -l, J ¥>(A|£„)AdA, 
and the diffusion coefficient as 

D(E,,t)=t^ J <fi(A\Ea)A^dA. 


(14) 


(15) 


Introducing the time-dependent variable E' =E + Ct reduces Equa¬ 
tion (13) to Einstein’s well-known diffusion equation 


dp 

IJi 



A=0 


(16) 


By construction, the initial conditions of the ensemble demand 
that p(E,t —>■ t(f\Ea,to) = S(E-Ea). This problem coincides with 
Einstein’s diffusion outwards from a point, where E-E^ plays the 
role of a ‘displacement’ from the initial energy E = Ea. The general 
solution to Equation (16) is a Green function (e.g. Krapivsky et al. 
2010) 


p{Ed\Ea,to)= , . = exp 

\/Alt D{Ea,t) 


[E-E,+C(Ead)f \ 

Ab(Ea,t) J’ 


with C = Ct' and D = Dr'. The mean and the dispersion of the dis¬ 
tribution are A = E-Ea= -C, and A^ - A = 2D, respectively. If the 
potential evolves adiabatically both coefficients C and D approach 
zero, and the probability function (17) becomes sharply peaked 
about E = Ea- In the limit R = R ^ 0, one has p(E,l\Ea,tQ) —> 
S(E—Ea), which recovers the adiabatic solution exactly. Section 2.3 
discusses the physical meaning of these coefficients and explores 
cases where C and D can be derived analytically. 

The evolution of a microcanonical ensemble in a time- 
dependent gravitational potential is therefore governed by Ein¬ 
stein’s equation for freely diffusing particles. A few aspects of this 
result are worth highlighting. First, notice that the distribution (17) 
is centred at A = E — Ea = —C. We shall see below that C for 
particle ensembles that have not in reached a state of dynamical 
equilibrium. Second, in contrast to the Brownian motion here the 
energy change A does not arise from a stochastic mechanical ex¬ 
change of energy with a granular medium but from the smooth time 
dependence of the gravitational field. This difference, albeit subtle, 
has very important consequences. Whereas in Einstein’s original 
paper one finds that the mean squared displacement of Brownian 
particles grows as = IDr' oc (TIK)t', where T corresponds to 
the temperature of the system and k to the viscosity coefficient, here 
the mean squared energy change with respect to the adiabatic solu¬ 
tion, A^, does not depend explicitly on r'. Indeed, the fact that r' 
vanishes from Equation (17) can be traced back to the continuous 
time-dependence of the gravitational field, which removes the prob¬ 
lem encountered by Einstein to describe the discontinuous nature of 
the Brownian phenomenon during very short time intervals^. The 
distinction between stochastic and deterministic processes needs to 
be emphasized, as the above results indicate that the evolution in 
phase space of systems subject to long-range forces will depend 
on whether the energy change is driven by a continuous dynamical 
process or by a random stochastic one, even though the equations 
that describe the evolution share the same form, i.e. that of Ein¬ 
stein’s diffusion equation. 


^ Indeed, in Einstein (1905) the coefficient D oc {x^)It' diverges in the 
limit r' —>■ 0. It became clear later on that during very short time intervals 
the Brownian motion does not behave like a Markovian process. 


2.3 Drift and diffusion coefficients 

Let us now turn to the thorny task of calculating the coefficients 
appearing in Equation (17). An analytical derivation of C and D 
is generally not possible due to the coupling between Equation (3) 
and the equations of motion (1), which introduces a dependence 
in the scale factor on the initial conditions of the orbit. However, 
for isotropic particle ensembles orbiting in scale-free gravitational 
fields it is possible to derive approximate solutions that offer useful 
insight into the physical meaning of these coefficients. In particular, 
it will be shown below that both C and D are closely related to 
quantities appearing in the virial theorem. 

We start by recalling that the probability function ip(A,Ea) 
is proportional to the number of particles in phase-space experi¬ 
encing an energy change E-Ea = A(r,v,r). Hence, Equations (14) 
and (15) can be written, respectively, as statistical averages of A 
and A^ over the phase-space volume ljj{E, t) accessible to particles 
with energy E, as discussed in Appendix A. Using Equation (A3) 
the coefficients C and D can be calculated as 

C(£,0«-- f iR/R)(r-y)S(E-H)d\d\, (18) 

w J 

and 

b{E,t) (R/Rf(r-yf5(E-H)d\d\, (19) 

where we have neglected all terms at order R and R^ in Equation (5), 
so that A = E — Ea ~ (R/R){r ■ v). This approximation is fairly ac¬ 
curate in potentials that do not evolve rapidly with time (see the 
numerical tests of P13). 


2.3.1 Scale free potentials 


In scale-free gravitational fields the scale factor R is independent of 
the phase-space trajectory of individual particles and only depends 
on the time f, hence Equations (18) and (19) can be simply written 
as C ~ -(R/R){r ■ v), and D ~ {R/Rf{{r ■ vf), with the brackets 
denoting microcanonical phase-space averages (see Appendix A). 
Note that C is proportional to (r • v), a quantity that is typically 
known as the virial function. This calls for the definition of a spe¬ 
cific moment of inertia 


T(E,t) 


1 

luj 


r^5(E-H)d?rd^\ 



r\2{E-A>)f'^dr, 


and a mean kinetic energy per degree of freedom 


( 20 ) 


T(E,t) 


1 f 

- / —S{E-H)d\d\ 
w / 2 


(4^ 

6uj 



r\2fE-^)f'^dr, 


( 21 ) 


where rm(E,t) is the maximum radius that particles with energy E 
can reach, that is T>(r„,r) = E. 

From the definition of the specific moment of inertia it is 
straightforward to show that T = {P/2) = {v ■ y) and T={v^) + 

(r • F) = 27’ -H W, where W{E,t) is the averaged potential energy of 
the microcanonical distribution (e.g. Heggie Sc Hut 2003). A state 
of dynamical equilibrium is therefore attained when distribution in 
configuration space obeys the virial theorem, T = T = 2T + W = 0. 

The drift coefficient (18) can be written in terms of virial quan- 
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n < -1 


n = -\ 


n> —1 


Phase-space volume p.u. energy (w) 
Constant (Bn) 

Specific moment of inertia (X) 
Mean kinetic energy (T) 


-(n+1) r[l/2-3/(n+l)] i ^ 1 1 /2„3 
3 7+n r[-3/(n+l)] 1^1 ' fn 

2 («+l)(7+n) r[l/2-3/(n+l)] r[2-3/(«+l)] 

3 7-3« r[3/2-3/(«+l)] r[-3/(«+l)] 

7+R r[l/2-5/(n+l)] r[-3/(n+l)] 2 

10 ii+« r[i/2-3/(«+i)] r[-5/(«+i)]'"f 
3 r[3/2-3/(«+l)] r[-3/(At+l)] I pi 
5-n r[l/2-3/(«+l)] r[2-3/(«+l)] 1^1 


(47r)2 



1/2 3 

' m 


1 

2 


(f 


3 ^ 3/2 2 
5/ 


a 

6 


r[l+3/(n+l)] |r|l/2..3 
6 r[3/2+3/(«+l)] 1^1 'fn 

2r[5/2+3/(o+l)] r[3/2+5/(n+l)] 
r[3/2+3/(«+l)] r[5/2+5/(R+l)] 

1 r[l+5/(»+l)] ri3/2+3/(n+l)] 2 

2 r[l+3/(«+l)] r[3/2+3/(R+l)] 

1 r[3/2+3/(»+l)] |p.| 

2 r[5/2+3/(«+l)] 1^1 


Table 1. Quantities defining the isotropic coefficient D{E^t) in Equation (24) derived in a gravitational field F(r,r) = Here E is the energy of the 

microcanonical ensemble, and rm{E,t) is the maximum radius that the particles with this energy can reach, that is = E, with ^ given by Equation (7). 

The Gamma function is defined as r(z) = exp[-?]dt. 


tides as 

C{E,t)=-(R/R)i, (22) 

which shows that the diffusion process outlined in §2.2 is closely 
related to the virializiation of systems out of equilibrium. Indeed, 
taking the time derivative on both sides of Equation (22) and ne¬ 
glecting terms at order R and R^ yields 

dC .... 

— « -(R/R)I = -(R/R)(2T + W). (23) 

at 

Microcanonical ensembles in virial equilibrium at f = to obey T ~ 
I ~ 0, which is equivalent to the condition C ~ dC/dt ~ 0. Taylor 
expanding the coefficient C at the time t = to+ t' therefore yields 
C(f) « C(to) + {dC/dt)T' « 0, which substituted in Equation (17) 
leads to E(t) « Ea{t). This is an important result, as we find that 
systems that are initially virialized follow an evolution that is close 
to the adiabatic limit at all times. We shall say that these systems 
evolve in quasi-dynamical equilibrium. 

Unfortunately, the derivation of C in systems that are neither 
phase mixed^ nor virialized at t = fo requires precise information 
on the trajectories of individual particles in phase space, which in 
general renders the problem analytically intractable. Tidal streams 
provide a particularly handy example of systems that are not in 
dynamical equilibrium and do not mix efficiently (see discussion 
in P13). Let us further illustrate the relation between phase mix¬ 
ing and virialization by discussing different states of a system at 
t = tQ (we shall come back to these examples in §4). First, one can 
distribute particles in phase space such that the number of them 
moving away from the centre of the potential equals the number 
moving toward it. By construction, this ensemble is in dynamical 
balance, i.e. (r • v) = T = 0 at t = to. Whether or not the system is 
also in virial equilibrium depends on the ratio between the kinetic 
and potential energies. As discussed above, if the system is initially 
virialized (I = 2T+W = 0) the evolution occurs in the limit of quasi- 
dynamical equilibrium, so that C ~ 0 at 1 = 1q + t' . However, if the 
condition for virial equilibrium is not satisfied (i.e. I ylO) the sys¬ 
tem will evolve away from its initial configuration towards an equi¬ 
librium state. During the process the drift coefficient C will exhibit 
periodic fluctuations which will damp progressively as the system 
regains its dynamical balance. A very similar process will be ob¬ 
served in systems whose initial configuration is not phase mixed 
(for example in tidal streams, where all particles are released from 

^ In this paper the term ‘mixed’ is used in a fairly broad context. In partic¬ 
ular, we shall say that particle ensembles are unmixed if they do not sample 
the accessible phase-space volume. An equivalent statement is that mixing 
is complete when phase-space occupation is maximum. Note that (i) the 
process of mixing happens in static as well as in time-dependent potentials, 
(ii) unmixed ensembles in general do not obey the ergodic principle. 


a progenitor galaxy with a wide range of velocities). In both cases 
collisionless relaxation proceeds in such a way that lim,_i.oo C = 0. 

In particle ensembles that are phase mixed the diffusion coeffi¬ 
cient D can be analytically expressed as function of virial quantities 

b(E,t) = Bn{R/RfTT, (24) 

where 6„ is a positive constant that only depends on the power- 
law index of the force (6). Table 1 provides analytical expres¬ 
sions for Bn, T and T for different power-law indices. The coef¬ 
ficient D is always positive and has a well-defined radial depen¬ 
dence which can be straightforwardly derived by inserting the val¬ 
ues of Table 1 into Equation (24). The kinematic energy scales as 
T ocE = <I>(r„) oc and the specific moment of inertia as T oc r^, 
hence the diffusion coefficient has a well-defined radial dependence 
b oc r^+". For indices n > -3, which covers the range of poten¬ 
tials of astronomical interest, b becomes very small at the centre of 
the potential and increases monotonically away from it, which im¬ 
plies that orbital diffusion predominantly affects the outer regions 
of gravitating systems. 


3 TIME-EVOLUTION OE THE DISTRIBUTION 
EUNCTION 

Section 2.2 relies on the construction of energy invariants to derive 
the evolution of a microcanonical ensemble of particles. Here we 
use those results to describe the non-equilibrium state of a system 
composed of an infinite number of microcanonical subsystems. 

Consider a system whose initial state is defined by the dis¬ 
tribution function f(Eo,{a}v,to), where {a}^ = ai,...,a^ are 
t 2 -integrals of motion related to the symmetry of the potential. 
If the normalization is chosen such that f /d^rd^v = 1, then 
/(EoirvljOiir, v],...,Q;iy[r,v],to)d^rd^v defines the probability to 
find a particle in the phase-space volume d^rd^ v centred at the coor¬ 
dinates (r, v) at the time to- In what follows it is convenient to work 
in the space of the integrals of motion. To this end we introduce the 
probability function 

N{E,{a}„,t) = f(E,{a}^,t)uj(E,{a}^,t), 

where uj is the density of states (see Appendix A for details). 

Our goal is to derive N(E, {a}v,t) at t = to-Hr' from the initial 
distribution AfEo,{a}i^,fo) in a time-dependent potential “Tfr,?). 
This shall be done in two steps. In Section 3.1 we derive an invari¬ 
ant probability distribution from the initial distribution function. In 
Section 3.2 we use the invariant distribution in order to find the dis¬ 
tribution of particles in the original integral-of-motion space at an 
arbitrary time 1 = to+ t'. 
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3.1 Invariant Distribution 


tion of Na requires one extra Green convolution 


In Section 2.1 we highlighted the fact that the canonical transfor¬ 
mation that removes the explicit time-dependence from the equa¬ 
tions of motion leaves the angular momentum invariant. This result 
implies that in time-dependent potentials which preserve the initial 
spatial symmetry the space defined by the constants of motion has 
u + l dimensions, namely the invariant energy I plus the t/-integrals 
a,. The Jeans theorem states that the probability function expressed 
as a function of the constants of the motion {/, ai,..., yields a 
steady-state solution of the collisionless Boltzmann equation, i.e. 
dA^ /At = (dN/dl)l + /daiy)di = dN/dt = 0. An equivalent 

statement is that the probability function N(I, tvi,..., a^,fo) remains 
invariant under the transformation to i—f = foH-r'. In what follows 
we emphasize the latter property by removing the explicit time- 
dependence from the invariant distribution, which shall be denoted 
as Al(/, {ajiy). In this work we are particularly interested in non¬ 
rotating spherical systems, which sets the integrals ( 01 , 02 , 03 ) = 
(Lx,Ly,L/), noting in passing that in potentials with time-varying 
spatial symmetry the diffusion process outlined in §2.2 takes place 
in multiple dimensions (see discussion in §5). 

Let us now derive the invariant distribution N(I,{o}v) from 
the initial probability function N(Eo, {a}^,ro). Using Equation (5) 
the invariant energy can be written as / = Rq(Eq - A), where Rq = 
R(tQ). Following the same steps as in Section 2 it is straightforward 
to show that the probability that a particle with an energy Eq at the 
time t = to has an energy in the interval /,/-td/ is 


p(I\EQ,to) = 


1 


47r£)(£o,to)J^o 


: exp 


[I-EoRl-C(Eo,to)Rlf 

Ab(Eo,to)R/, 


(25) 

with the coefficients C and D given by Equations (14) and (15), 
respectively. 

The probability function p(l\Eo,to) is a solution to the differ¬ 
ential Equation (16) and therefore a Green function (e.g. Krapivsky 
et al. 2010). Hence, the invariant probability function is found by 
convolving the initial distribution function N(Eo,{a},j,to) with the 
probability function (25) 


A'(/,{a}.) = j p(I\Eo,to)N(Eo,{a},x,to)AEo (26) 

_ r N(Eo,{o},,to) / [I-EoRl-C(Eo,to)Rlf\.^ 

^4nD(Eo,to)R‘*o ^ 4D(Eo,to)Ro / 

with the limits of the integral defined by the range of energies of 
the system (see examples in §4). 


3.2 Time-dependent Distribution Function 

One can now derive N(E, {a}v,t) from the initial probability func¬ 
tion N(Eo, {o}v,to) following similar steps as as in §3.1. 

Recall that the energy can be written as £ = Ea(t) + A, where 
Ea = I/iP is the adiabatic energy. Let us introduce the function 
Na(Ea,{o}v,t) which defines the probability to find a particle 
within the adiabatic energy interval Ea,Ea + AEa at the time r. In po¬ 
tentials with a fixed symmetry the probability function at the time 
t = to + t' will be the result of the following Green convolution 

N(E,{a},x,l) = j p(E,t\Ea,to)Na(Ea,{o},x,to)AEa, (27) 

where the Green function p is given by Equation (17). The deriva- 


^a(Ea, {fV }[/ , t) 


p(E„t\I)N(I,{o},x)dI. 


(28) 


In scale-free fields the probability distribution p(Ea,l\I) = 5(Ea — 
I/lP) because R does not have an intrinsic orbital dependence (see 
§2.1.1). In this case Na can be directly obtained from the invariant 
function (26) as 


Na(Ea,{a},x,t)AEa=N(I[Eal{o},x) 


AEa 


R^AEa=N(I,{o},)AI. (29) 


Substituting (29) into (27), inserting Equation (26) and re-arranging 
the integrands yields 


iV(£,{a}.,t) = J p(E,t\IR-\to)N(I,{a},x)AI 
AEoN(EQ,{a},x,lo) J p(E,t\IRr^ ,to)p(I\Eo,to)dI 
= J AEoN(Eo,{o},,,to)pciE,t\EQ,to). 


(30) 


The probability function pc(E,t\Eo,lo) is the result of the convolu¬ 
tion of two Gaussians. In general, this function cannot be expressed 
analytically owing to the non-trivial dependence of the function 
p(E,t\IRr^,to) on the energy /. However, in potentials that evolve 
slowly the function pc will peak sharply around E ~ Ea = I /PJ. Re¬ 
placing the coefficients C(IRr^,t) and b(IRr^,t) on the right-hand 
integral of Equation (30) by C(E,t) a.nAD(E,t), respectively, yields 
an analytical expression for pc which, not surprisingly, also is a 
Gaussian function 

Pc(E,t\Eo,to) = J p(E,t\IRr^,to)p(I\EQ,to)AI (31) 

1 r [E-Eo(Ro/R)^ + ACf 

~ , . exp-^-=- 

V 4ttDc I 

where 

AC(E,t,Eo,to) = C(E,t)-C(Eo,to){^'^ , (32) 

and 

ba(E,t,Eo,to) = b(E,t) + b(Eo,to){^'^ . (33) 

Section 4 shows that non-equilibrium systems wherein AC ^ 0 ex¬ 
hibit periodic fluctuations in the distribution function which damp 
out progressively as the system approaches an equilibrium state. 

Notice that pc plays the role of a transition probability be¬ 
tween the states to and to + r'. Appendix B discusses how transi¬ 
tion probabilities can be used to describe the evolutionary path be¬ 
tween time intervals of arbitrary length through the construction of 
Markov chains. 


3.3 The Fokker-Planck approximation 

This Section relates N(E,{o},x,t) to the distribution obtained in 
the adiabatic limit, Na(Ea,{o}v,t) by finding an approximate so¬ 
lution to Equation (27). To this end we apply a technique that is 
commonly used in variational calculus, where one assumes that the 
transformation Na 1 -^ N takes place during a short, but non-zero, 
time interval to,lo + T'. The relation between the two functions at 
an arbitraty time to is made by Taylor expanding N and taking the 
limit t' —>■ 0. 

Recall that in a slowly-varying gravitational potential the 
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Green function p(£,f|£'a,ro)willbe sharply peaked around E =Ea, 
which calls for using A = E—Ea as the integration variable in Equa¬ 
tion (27). Since A is given as a function of phase-space coordi¬ 
nates one needs to replace the probability function p{E,t\Ea,t(i) 
by ip(A\Ea). In addition, given that {a}^ are assumed to be v- 
constants of the motion we can simplify our notation by marginaliz¬ 
ing over J d'^a on both sides of Equation (27), which then becomes 

N{E,to+r') = [ p(A\E-A)Na(E-A,tQ)dA. (34) 


For an arbitrarily small r' one can now follow the exact same steps 
as in §2.2. We expand the left-side of Equation (34) as a Taylor 
series in r' up to the first order, and similarly expand the function 
ipN within the right-hand integral as a Taylor series in A to the 
second order, which yields 


Na(E,to)+T'^+... = j dA[y>-Vt/sA-t ^vVA^ + ---] 
x[iV.-ViV„A-H^V'iV„AV...] 


N 

^ dA[ipNa-V{pNa) A 

+ ^V^(ipNa)A^ 

! 99 dA-V 

Na J pAdA 

+ 

2 

Na J pA^dA 


(35) 


where we have defined the operator V = d/dE and used the initial 
conditions N(E,to) = Na(E,to). 

From Equations (9), (14) and (15) one has that f ipAA = 1, 
f (pAdA = -C and f ipA^dA = 2D, respectively. Hence, after elim¬ 
inating Na{E,to) from both sides of Equation (35) we find 


r' — « V(CN,) + V^(bN,), (36) 

at 

which corresponds to the Fokker-Planck equation formulated in the 
energy space (e.g. Spitzer 1987). 

Given that the right-hand side of Equation (36) is independent 
of the value of r' (see discussion in §2.2) the time interval can be 
made arbitrarily small. By taking the limit r' —>■ 0 we can write the 
time derivative as 

dN _ N(E,to)-Na{E,to) 
dt ~ t' 

which reduces Equation (36) to 


Al ~ H- V(C%) -H V^{bNa). 


(37) 


Equation (37) therefore provides an approximate solution to (27) 
after marginalizing over d'' a. A few points of interest are worth dis¬ 
cussing. First, in potentials that evolve adiabatically the coefficients 
become C = Zf ~ 0 (see §2.2), which leads to A ~ TVa. This recov¬ 
ers the well-known fact that in slowly-varying systems the shape 
of distribution function f = N/w is an adiabatic invariant. Second, 
notice that the first- and second-order corrections beyond the adia¬ 
batic limit describe a net flux of particles in the energy space. From 
Equation (36) the probability flux crossing the energy layer E can 
be estimated as 

J(E) = j _ cNa + \/(DNa). (38) 

Not surprisingly, the flux diverges in the limit r' —> 0 because in 
that limit the energy change becomes instantaneous. However, the 
bulk displacement of particles in energy space, J = Jt = CNa + 
\7(DNa), is again independent of the choice of r'. 

In scale-free fields it is possible to derive a straightfoward 
interpretation of Equation (38). At first order the flux scales as 


J ~ CNa oc (R/R){r ■ \)Na, with the brackets denoting average over 
microcanonical ensembles with energy E (see Appendix A). The 
flux, therefore, is the result of internal macroscopic motions and 
oscillates in phase with the net radial velocity of the particle en¬ 
semble. If the system is in a state of quasi-dynamical equilibrium 
then (r • v) « 0, and therefore C « 0. In this limit the flux becomes 
diffusive, /« V(DAj). The fact that Z) > 0, VZ7 > 0 and \E\ oc T 
(see Table 1) implies that particles will diffuse in a direction oppo¬ 
site to the temperature gradient. In gravitationally-bound systems 
this generally leads to a net particle flux from regions of high den¬ 
sity/temperature toward regions of low density/temperature. 


4 A-BODY EXPERIMENTS 


In this Section we run a number of restricted A-body experiments 
that illustrate the evolution of tracer particles orbiting in the time- 
dependent scale-free potential given by Equation (7). To simplify 
the analysis we shall consider spherical power-law potentials that 
evolve linearly with time as 


Uh) _ , t-tp 

“ ^ ^ rt ? 

Mo -^0 


(39) 


where Pp = P(E, tp) = 2 drj^2(E - <I>[r, fo]) is the radial period 

of an orbit with /■„ = /io = 1 at r = to- 

Particles are distributed within a distance range from the cen¬ 
tre of the potential, rmin < r < rpm, with rmin = 0.1 and rum =10. The 
initial radii and velocities of the particle ensemble are drawn from 
an initial distribution function f(Ep,tp) using a rejection algorithm, 
while the position and velocity vectors are isotropically distributed 
over a sphere. The equations of motion (1) of each individual parti¬ 
cle are then integrated forward in time using a leap-frog algorithm 
whose time-step is chosen such that the energy is conserved at least 
at a lO^"* level in the static case (e = 0). 


4.1 Analytical diffusion coefficients 


First we shall study a few examples where the coefficients C and 
D can be derived analytically using the dynamical invariants intro¬ 
duced in §2.1. 

Let fiEp.tp) be a distribution function where the number of 
particles in the energy interval E,E + dE has a fixed value indepen¬ 
dently of E. From Appendix A this condition implies 


f{Ep,tp) 


Np 

uj{Ep,tpy 


(40) 


Integrating on both sides of Equation (40) over the phase-space vol¬ 
ume d^G shows that the total number of particles in the model is 
NpAE, where AE corresponds to the range of energies covered by 
the particle ensemble. 

Fig. 1 shows the variation of energy relative to the adiabatic 
value, i.e. E—Ea{t), as a function of the orbital energy in time- 
dependent potentials with a growth rate e = -0.1. The left, mid¬ 
dle and right panels adopt Keplerian (n = -2), logarithmic (« = -1) 
and harmonic (n = +i) potentials, respectively. In general, we find 
that the displacements from the adiabatic limit are on average very 
small, E—Ea~Q, which indicates that the A-body models are ini¬ 
tially in dynamical equilibrium. In contrast, the variance (E-Ea)^ 
increases significantly as one moves away from the centre of the 
potential (i.e. toward the right-hand side of each panel). This be¬ 
haviour can be easily understood using the diffusion coefficients 
derived in §2.3.1. In a scale-free gravitational field (E —Eay = 2D oc 



8 Jorge Penarrubia 



-3 - 2 - 10 - 3 - 2-10 1 2 3 


E/Mo 



E/ Mo 


Figure 1. Upper panel: Energy change with respect to the adiabatic value of N* = 10"^ particles orbiting in a scale-free gravitational potential (7) that grows 
linearly as p,{t)/pQ = 1 +e{t — tQ)/PQ, where /io = p{to) and Pq is the radial period of an orbit with energy = O.lrjim at r = ?o (see text). These models are 
run from t = to to t = to+ Pq with e = —0.1. All models show a mean energy C = E—Ea ^ 0, which indicates that their evolution proceeds in quasi-dynamical 
equilibrium. In these three cases the energy dispersion increases as one moves away from the centre of the potential (i.e. toward the right-hand side of each 
panel). Lower panel. Energy variance as a function of energy for particles orbiting in the potentials of the upper panels. Different symbols correspond to 
potentials evolving at different rates. Note that the formula {E—Ea)^ = + 2D ^ 2D, with D derived from Equation (24) and Table 1 (dashed lines), is in 

excellent agreement with the numerical results. 


where rm(E^t) is the maximum radius that a particle with an 
energy E can reach in the potential (7), that is 


rm{E,t) 


nimexp[£/Mf)] = 


(41) 


The constant in front of the potential (7) is <l>oo = 0 for the Keplerian 
and harmonic cases, and $oo = -'l’(riim,ro) for the logarithmic one. 
Hence, from Equation (41) it follows that D oc £ ', D oc and 

D xE^ for power-law indices n = -2, -1, and -Hi, respectively. The 
lower panels of Fig. 1 shows that the diffusion coefficients given 
by Equation (24) and listed in Table 1 (dashed lines) provide an 
excellent match against the A^-body models (symbols) over a large 
range of energies and growth rates. 


4.2 Violent relaxation 

In this Section we study the evolution of A^-body models with a 
virial ratio |2r /W\ <C 1 at f = fo. Setting the initial kinetic pres¬ 
sure to such a low value leads to a rapid infall of particles toward 


the centre of the potential, a process typically known as cold col¬ 
lapse. Given that the potential varies as a function of time the or¬ 
bital energies of individual particles change in a non-trivial way 
as the system approaches equilibrium (Lynden-Bell 1967; Mo, van 
den Bosch & White 2010). 

As an illustration of this process we generate A-body equilib¬ 
rium realizations of a lowered Maxwellian distribution in a loga¬ 
rithmic potential 


/(£o,fo) 


A\e -l] ,£o < l&iim, 

0 , £o ^ (hlim 5 


(42) 


where A is a normalization factor, and $1101 = po In(nim) is an energy 
truncation. To simplify our models we set = po/2 = 3T, where 
T is the kinetic energy associated with the potential (see Table 1). 

Next, we multiply the velocities of all particles by a factor 
q. To highlight non-equilibrium features we choose a small value, 
q = 0.2, which leads to |2r /W\ = q^ = 0.04. Such a low virial ratio 
guarantees the collapse of the system on a time-scale comparable 
to its free-fall time. Fig. 2 shows ten snap-shots of a model orbiting 
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Figure 2. Cold collapse of a model with |27'/VF| = 0.04 at / = fQ. Particles move in a time-dependent potential that vaiies at a constant rate e = —0.1. The system 
as a whole evolves toward an equilibrium configuration following a ‘violent relaxation’ process. 



Figure 3. Drift coefficient C = -{R/R){r-\) as a function of energy of the model shown in Fig. 2. Notice that at f = to the system is phased mixed (i.e. (r ■ v) = 0) 
but out of virial equilibrium, |2r/VF| = 0.04. By the end of the simulation only the internal regions of the potential have reached a state of quasi-dynamical 
equilibrium (C ~ 0). 


in a time-dependent potential that evolves at a rate e = -0.1. Cold 
collapse happens early on (r' « 0.2Po) and leads to the formation 
of shell structures that move progressively towards larger radii with 
time. By the end of the simulation, t' ~ 2.0Po, the model has not 
yet reached dynamical equilibrium. By definition the system is in a 
state of ‘incomplete relaxation’'*. 

Fig. 3 shows the drift coefficient C = -(R/R){r ■ v) as a func¬ 
tion of energy for the snap-shots shown in Fig. 2. The initial Al-body 
models are phase mixed, such that C oc (r ■ v) « 0 for all energies. 

These models provide a useful representation of the dynamical state of 
the outer regions of galactic haloes and galaxy clusters, where dynamical 
times are comparable to the age of the Universe and phase mixing becomes 
very inefficient (e.g. Henon 1964, Lynden-Bell 1967; Dehnen 2005). 


Flowever, the fact that \2T /VF| 1 implies that the model is far 
from a virialized state. As the system begins to collapse the aver¬ 
aged radial velocity of the particle ensemble is negative at all radii 
and dC/dt > 0 at all energies. At slightly later times, r' > 0.2Po, 
the coefficient C begins to exhibit coherent fluctuations in the inner¬ 
most regions of the potential (left side of the panels), as particles 
with short orbital period go through pericentre and start moving to¬ 
ward larger radii. In the outskirts, however, particles are still falling 
in from large distances, which translates into positive values of C. 
The negative crests are associated with the shell features of Fig. 2. 
Given that in a potential with n = —\ the orbital period decreases 
toward the central regions of the potential as 

P{Eq) = (27r//r)‘^^?-„, = (27r/^)‘''^f-iimexp(£//i), 
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Figure 4. Upper panels: Distribution function / = NJu) of the models shown in Fig. 2 at three different snap-shots. Red solid lines correspond to the Green 
convolution given by Equation (30) with coefficients C(E,t) = —(R/R){r ■ v) and D(E,t) = {R/R)^{{r ■ v)^) measured from the phase-space locations of the N- 
body particles. As expected, the distribution /(/) of Equation (26) (cyan dashed curves) remains invariant throughout the evolution of the system. Lower panels: 
Deviation of the A-body models from the adiabatic distribution, AA(, = N(E,t)-Na{E,t) (green dots). The Eokker-Planck approximation. Equation (37), is 
shown with magenta solid lines. Note that the ripples and troughs of the curves are located at energies where the gradient V(CA) finds local maxima and 
minima, respectively. Departures from the adiabatic solution are particularly strong in the outskirts of the system, E > 


fluctuations in C damp out progressively from inside out. By r' = 
1.9/3) the mean radial velocity of particles with E — < -2 is 

(r • v) « 0, signalling that the inner regions of the system are phase- 
mixed and evolving in a state of quasi-dynamical equilibrium. 

Fig. 4 illustrates the complexity involved in describing the 
state of systems undergoing violent relaxation. Green dots show the 
distribution function of the A^-body models plotted in Fig. 2 at three 
different snap-shots. The decreasing potential (e < 0) shifts the or¬ 
bital energies to higher values, which leads to a non-monotonic 
increase of f(E,t) at fixed energies. Interestingly, the distribution 
function is not completely smooth. Non-equilibrium features arise 
in the inner-most regions of the potential and propagate toward high 
energies as time goes by. After constructing a larger suite of A-body 
models (not shown here) we find that the amplitude of these fluctu¬ 
ations increases for larger e and smaller q values, which correspond 
to faster growth rates and higher radial anisotropies, respectively. 

Equation (30) (red solid lines) is able to capture these com¬ 
plexities and provide an accurate statistical description of the non¬ 
equilibrium state of the system. For simplicity we assume that the 
transition probability pc{E,t\Eo,to) independent of the evolu¬ 
tionary path taken by the system between ?o and ro + t' (see Ap¬ 
pendix B). The Green convolution is solved by setting C(Eo,to] = 0 
and computing D(Eo,tQ) analytically from Equation (24) and Ta¬ 
ble 1. The coefficients C(E.t) and D(E,t) are measured from the 
phase-space coordinates of the A-body particles at ? = fo H- r' as 
C(E,t) = -{R/R){r-\} and D{E,t) = {R/Rf {(r ■ \f}. Note that in 
this case the solutions to Equation (30) are not purely statistical, as 
our calculation of the coefficients adds additional dynamical infor¬ 
mation. 

The Eokker-Planck approximation provides useful insight into 


the dynamical mechanisms that drive the shape of the distribution 
function. The lower panels of Eig. 4 plots the difference between 
the A-body distribution function and the adiabatic solution given 
by Equation (29) at fixed energy values (green dots). From Equa¬ 
tion (39) particles respond adiabatically to a time-varying force 
if their orbital periods obey P(tE>/<l>o) = e(P/Po) <C 1. Since P in¬ 
creases exponentially with the particle energy we find the distribu¬ 
tion A evolves adiabatically (AAa ~ 0) at /? ^ 4>um, while strong 
departures from the adiabatic solution (| AAa| ~ A) are visible at 
E > tDiim, where the potential changes significantly during an or¬ 
bital period, i.e. P($/$o) ~ 1. As a result, the Eokker-Planck ap¬ 
proximation (solid magenta lines) becomes less accurate in the out¬ 
skirts of the system. Note that similar deviations are also visible in 
the upper panels at r' = 1.8Po, albeit with a lesser magnitude. 


The existence of internal macroscopic motions (C ^ 0) leads 
to fluctuations of the distribution function which travel toward high 
energies, as shown in Pig. 3. Comparison with the solid magenta 
lines shows that the ripples and troughs of the fluctuations are lo¬ 
cated at energies where the gradient V(CA) = d(CN)ldE finds lo¬ 
cal maxima and minima, respectively. Notice that by the end of the 
simulation (r' = 1.8Po) relaxation is still ‘incomplete’. 


In phase-mixed regions relaxation proceeds in quasi- 
dynamical equilibrium (C ~ 0), and the evolution of the distribution 
function becomes diffusive, as discussed in §3.3. Tho final state that 
emerges from violent relaxation can be found by setting C{E, t) = 0 
in Equation (30). 












Non-equilibrium gravitational mechanics 11 


5 SUMMARY & DISCUSSION 

This paper introduces a probability theory that describes the non¬ 
equilibrium evolution of large particle ensembles orbiting in a time- 
dependent gravitational field. The theory is constructed on the ba¬ 
sis of dynamical invariants (quantities conserved along the phase- 
space path of a particle motion) and does not rely on maximum- 
entropy probabilities or ergodicity assumptions. 

A fundamental result of this paper states that the evolution of 
particles with the same orbital energy (the so-called microcanoni- 
cal ensemble) has the same form as Einstein’s equation for freely 
diffusing particles with coefficients that are closely related to virial 
quantities, such as the moment of inertia and the temperature of the 
system. Since the general solution to the diffusion equation corre¬ 
sponds to a Green function, the non-equilibrium distribution func¬ 
tion of systems composed of an infinite number of microcanoni- 
cal subsystems is found by convolving the initial distribution func¬ 
tion with the solution to the diffusion equation. The Fokker-Planck 
equations (37) provide an approximate solution to this type of 
Green convolutions. Thus, we conclude that the problem of col¬ 
lisionless relaxation in the linear regime reduces to the derivation 
of diffusion coefficients. 

The analogy between our mathematical formalism and 
stochastic calculus is undeniable and also somewhat striking. In¬ 
deed, at first glance it is not completely obvious why equations 
that are common to a large range of stochastic processes do also 
describe systems governed by fully deterministic equations of mo¬ 
tion. To answer this question it is useful to examine the role of the 
probability function ifi(A) in some detail. In Einstein (1905) the 
explicit form of v?(A) contains all information on the dynamics of 
the collisions that drive the motion of Brownian particles. How¬ 
ever, its role is that of a nuisance function, for the diffusion co¬ 
efficients are derived indirectly from observations of macroscopic 
quantities (such as the mean squared displacement of the particles 
as well as the temperature and the viscosity of the medium where 
the particles are suspended) over an extended period of time, r'. 
The great virtue of this approach is that it permits a derivation of 
the probability distribution of random displacements at the time 
tQ = to + r' from the distribution at t = to without having to spec¬ 
ify the form of (p(A)! Although this frees us from the necessity 
to understand the mechanics of Brownian motion at a microscopic 
level, it becomes necessary to assume that time-averaged proper¬ 
ties emerge from phase-space averages of microscopic quantities, 
which constitutes the crux of the ergodic principle. Unfortunately, 
the equivalence between time average and average over ensembles 
only arises when the system can visit all the possible microstates, 
many times, during a long period of time, which in general can¬ 
not be guaranteed in non-equilibrium systems. In contrast, in the 
present work the form of (p(A) is known and corresponds to the 
microcanonical distribution function (A4), while the dynamical in¬ 
variants of Section 2.1 define the phase-space dependence of A. In 
systems with known dynamical invariants it is therefore possible 
to find self-consistent solutions to the diffusion equation without 
resorting to time averages of macroscopic quantities, thus remov¬ 
ing the dependence on the time interval r' from the solution and - 
more importantly - the necessity to rely on ergodicity assumptions 
in order to describe the evolution of gravitating systems. 

Alas, not all potentials admit analytical solutions to the diffu¬ 
sion equation. For example, in systems that are not initially viri- 
alized a derivation of the diffusion coefficients requires precise 
knowledge of the trajectories of individual particles in phase space, 
rendering the problem analytically intractable. Also, to this date 


an analytical derivation of energy invariants is only feasible in 
scale-free gravitational potentials (see P13). And yet the theory has 
many interesting applications that go beyond the analytical cases 
explored here. For example notice that in this paper the time evolu¬ 
tion of the potential has been intentionally detached from the spa¬ 
tial distribution of the particles under scrutiny. In essence, we have 
treated particles as mass-less tracers of an underlying gravitational 
field which is completely specified by the particle coordinates and 
the time (that is, the field is purely mechanical and not a statistical 
object). Although there is a large number of astrophysical objects 
that meet these conditions (stars in dwarf spheroidal galaxies, tidal 
streams in galactic haloes, and galaxies in galaxy clusters are a few 
prototypical examples of kinematic tracers in time-dependent po¬ 
tentials), an exciting task for the future may be to study a more 
constrained case where the initial state of the system is specified by 
/(r, v,fo) at f = to, and the evolution of particle ensemble is driven 
by its own self-gravity 

V^$(r,0 = TttG [ /(r,v,r)d^v, 


with f = N/oj given by Equation (26). A solution to this problem 
will provide a deep insight into the evolution and stability of self- 
gravitating collisionless systems. 

By considering time-dependent potentials with a fixed spa¬ 
tial symmetry we have confined the diffusion process outlined in 
§2.2 to the energy dimension. Yet, our formalism can be straight¬ 
forwardly extended to systems with a varying spatial symmetry. As 
an illustration, let us consider a system in which both, the angular 
momentum and the energy of individual particles are allowed to 
vary with time. The Green function appearing in Equation (26) has 
now one extra dimension, i.e. p{E,L,t\E',L',to), which defines the 
probability that a particle with an energy E' and angular momen¬ 
tum L! at the time t = to has an energy in the interval E,E + AE and 
an angular momentum in the interval L,L+AL at the time t = to+ t' . 
Armed with the new probability function one can follow the same 
steps as in §3.3 to derive the two-dimensional Fokker-Planck equa¬ 
tions 

, OA d — d — 1 — 

where the coefficients correspond to phase-space averages of AL = 
L-L', AL^ = {L-L'f and A(AL) = (E-E')(L-L'), as indicated 
by Equation (A7). The simplest case corresponds to an evolution¬ 
ary process in which the change of angular momentum is inde¬ 
pendent of the orbital energy of individual particles. The cross¬ 
coefficient becomes A(AL) = 0, which leads to a Green function 
p{E,L,t\E',L',to) with a relatively simple form 


1 

--- exp 

[(47r)2£)£)z,]'/2 


(A-hC)" (AL + Ctf 
AD ADl 


where AL = -Cl{L' ,t) and AL^ — AlJ' = IDtiL' ,t) play the same 
role as C and D in the energy space. In general, radial orbit instabil¬ 
ities (see Henon 1973, Huss et al. 1999, MacMillan et al. 2006) will 
lead to A(AL) / 0 even in spherical potentials. The importance of 
orbital inestabilities has been recently highlighted by Pontzen et al. 
(2015), who argue that collisionless systems evolve toward attractor 
states which are maximally stable when external perturbations are 
applied. However, it is important to emphasize that most galaxies 
never reach a fully-relaxed state within a Hubble time, which may 
be the reason for the observed mismatch between the distribution 
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function derived from maximum-entropy arguments and those ob¬ 
served in cosmological A^-body simulations (Pontzen & Governato 
2013). Given that our theory can describe the non-equilibrium state 
of self-gravitating systems, a derivation of the crossed-coefficients 
A(AL) will shed light on the way in which the distribution function 
actually evolves towards its maximally-stable limit. 

Finally, the fact that the evolution of the microcanonical distri¬ 
bution has the same form as Einstein’s equation for purely stochas¬ 
tic processes offers the tantalizing possibility to incorporate the ef¬ 
fects of random particle-particle collisions into our probability the¬ 
ory in a natural way. Such an extension may be useful, for example, 
to study a range of dynamical processes taking place in planetary 
systems, dense stellar objects and/or self-interacting dark matter 
(SIDM) haloes. 


6 ACKNOWLEDGEMENTS 

I would like to thank Justin Read and Frank van den Bosch for their 
insightful questions, which proved to be key to clarify essential as¬ 
pects of the theory. Also, many thanks to Matt Walker for pointing 
me to Donald’s paper on Dirac’s cosmology, and to the DwaIN net¬ 
work for tete-a-tete discussions. 


REFERENCES 

Antonov, V. A. 1961, Soviet Ast., 4, 859 
Arad, I., & Lynden-Bell, D. 2005, MNRAS, 361, 385 
Arad, I., & Johansson, P. H. 2005, MNRAS, 362, 252 
Campa, A., Dauxois, T., & Ruffo, S. 2009, Phys. Rep., 480, 57 
Chavanis, P. H., Sommeria, J., & Robert, R. 1996, ApJ, 471, 385 
Dehnen, W. 2005, MNRAS, 360, 892 
Dejonghe, H. 1987, ApJ, 320, 477 
Einstein, A. 1905, Annalen der Physik, 322, 549 
Heggie, D., & Hut, P. 2003, The Gravitational Million-Body Prob¬ 
lem: A Multidisciplinary Approach to Star Cluster Dynam¬ 
ics, hy Douglas Heggie and Piet Hut. Cambridge University 
Press, 2003, 372 pp., 

Henon, M. 1964, Annales d’Astrophysique, 27, 83 
Henon, M. 1973, A&A, 24, 229 
Jaynes, E. T. 1957, Physical Review, 106, 620 
Joyce, M., Marcos, B., & Sylos Labini, F. 2009, Journal of Statisti¬ 
cal Mechanics: Theory and Experiment, 4, 19 
Huss, A., Jain, B., & Steinmetz, M. 1999, ApJ, 517, 64 
Krapivsky, P. L., Redner, S., & Ben-Naim, E. 2010, A Kinetic View 
of Statistical Physics, by Pavel L. Krapivsky , Sidney Redner, 
Eli Ben-Naim, Cambridge, UK: Cambridge University Press, 
2010, 

Landau, L. D., & Lifshitz, E. M. 1980, Course of theoretical 
physics, Pergamon International Library of Science, Technol¬ 
ogy, Engineering and Social Studies, Oxford: Pergamon Press, 
19801c 1980, 3rd rev.and enlarg. ed.. 

Levin, Y., Pakter, R., & Rizzato, E. B. 2008, Phys. Rev. E, 78, 
021130 

Levin, Y, Pakter, R., Rizzato, F. B., Teles, T. N., & Benetti, F. P. C. 

2014, Phys. Rep., 535, 1 
Lynden-Bell, D. 1967, MNRAS, 136, 101 
Lynden-Bell, D., & Lynden-Bell, R. M. 1977, MNRAS, 181, 405 
Lynden-Bell, D. 1982, The Observatory, 102, 86 
Lynden-Bell, D. 1999, Physica A Statistical Mechanics and its Ap¬ 
plications, 263, 293 


MacMillan, J. D., Widrow, L. M., & Henriksen, R. N. 2006, ApJ, 
653, 43 

Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation 
and Evolution, by Houjun Mo , Frank van den Bosch , Simon 
White, Cambridge, UK: Cambridge University Press, 2010, 
Nakamura, T. K. 2000, ApJ, 531, 739 
Padmanabhan, T. 1990, Phys. Rep., 188, 285 
Penarrubia, J. 2013, MNRAS, 433, 2576 (P13) 

Pontzen, A., & Governato, F. 2013, MNRAS, 430, 121 
Pontzen, A., Read, J., Teyssier, R., et al. 2015, arXiv: 1502.07356 
Robert, R., & Sommeria, J. 1992, Physical Review Letters, 69, 
2776 

Spitzer, L. 1987, Princeton, NJ, Princeton University Press, 1987, 
191 p., 

Sridhar, S. 1987, Journal of Astrophysics and Astronomy, 8, 257 
Stachel, J. 1998, Einstein’s miraculous year : five papers that 
changed the face of physics / edited and introduced by John 
Stachel Princeton, N.J. : Princeton University Press 
Sylos Labini, R, Benhaiem, D., & Joyce, M. 2015, MNRAS, 449, 
4458 

Wax, N. 1954, New York: Dover Publication, 1954, edited by Wax, 
Nelson, 


APPENDIX A: STATISTICAL AVERAGING 

The large number of particles typically enclosed in astrophysical 
systems, together with the necessary incompleteness of observa¬ 
tional data, makes phase-space averaging one of the most useful 
tools in Statistical Mechanics. Here we discuss the use of a micro- 
canonical distribution to calculate phase-space averages of macro¬ 
scopic quantities. 

In order to express the distribution function introduced in §3 
as a function of phase-space coordinates one needs to map points 
in the integral-of-motion space onto the phase-space volume d^O = 
d^rd^v through the equation 

/(r,v,f)d'’H = /(£',{a}^,r) -- AEiPa (Al) 

d{E,au...,au) 

= N{E,{a}^,mEA'"a, 

where [w] = d{T,y)/d(E,a\,...,av) is the matrix of density of 
states. The Jacobian of this matrix, w(E, {a}^,t), defines the max¬ 
imum space-space volume that particles with a given combina¬ 
tion of integrals of motion can potentially sample. With this def¬ 
inition the probability of finding a particle within the energy in¬ 
terval {E,E + AE) can be found by marginalizing Equation (Al) 
over {a}ly, that is N{E,t) = f A)/?, {a}iy,f)d'^a, with the probabil¬ 
ity function normalized such that f N{E,t)AE = 1. 

A very useful thermodynamical concept corresponds to mi¬ 
crocanonical ensembles. Consider a gravitating system where all 
particles have a known orbital energy. Now isolate a subset of par¬ 
ticles with a specific energy, E. Even if the integration of the equa¬ 
tions of motion were feasible, that information alone would not 
be enough to determine the trajectories of those particles in phase 
space, for it is not possible to define the initial conditions of the 
problem. This means that in general we will have a complete lack 
of knowledge of their microscopic state. The only available infor¬ 
mation reduces to the condition imposed on their phase-space lo¬ 
cations, which must be confined to a hyper-surface in phase space, 
El{E), where E = H(r, v, t), and H = l/2v^ -H $ is the Hamiltonian 
of the system. 
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Averaging X over all of phase space of the system yields 


A key question is, how are those particles distributed on the 
single-energy hypersurface fl(£)? If they are homogeneously dis¬ 
tributed one can easily construct an ensemble of microcanonical 
states by randomly sampling fl(£) a very large number of times. 
By construction all microstates would be equally probable, and in¬ 
tuition tells us that the averaged properties of a microstate would 
be indistinguishable from the rest within small fluctuations. But 
what do we mean exactly with ‘equally probable’? Here we can 
either adopt a frequentist definition, where one measures the av¬ 
eraged properties of each microstate at different times in order to 
assign a certain probability, or a statistical one, where the assigned 
probability is based on the density of particles on Q,{E). If the sys¬ 
tem is ergodic both probabilities are the same within small fluc¬ 
tuations, which immediately leads to an equivalence between time 
average and average over microstates. However, one can easily find 
examples of gravitating systems where ergodicity is broken. Con¬ 
sider, for example, a Globular Cluster losing stars to the host’s tidal 
field. Initially these stars will distribute along tidal tails with a small 
spread in orbital energy, which can be neglected for the sake of ar¬ 
gument. Clearly, since tidal tails do not sample Q.{E) in a randomn 
fashion the ergodic hypothesis would lead to a wrong description 
of their macroscopic properties. However, if we wait long enough 
we would see the tidal tails spreading throughout larger segments 
of this hypersurface, a process that is typically known as phase mix¬ 
ing. Thus, the system becomes ergodic in the limit t —>■ oo, as the 
occupation of Q.{E) evolves toward a global maximum (see P13 for 
numerical illustrations of this process). Although one can formally 
construct ensembles of ‘unmixed’ microstates (see §2), a descrip¬ 
tion of their macroscopic properties will require dynamical infor¬ 
mation on the trajectories of particles in phase space, which in gen¬ 
eral renders the problem non analytical, as discussed in §4. In the 
remainder of the Appendix we focus on phase-mixed systems, for 
which microcanonical averages can be straightforwardly derived. 

For spherical & isotropic systems the density of states of a 
microcanonical ensemble of particles with an energy E reduces to 
uj = dQ/dE, which can be written as 

uj(E.t) = e(E-H)(fQ= f 5(E-H)d‘'n (A2) 

oE J 

= (47r)^ [ [2(E - $)] ‘^V^dr, 

Jo 

wher O(E-H) the Heaviside function, and r^iE, t) is the maximum 
radius that particles with energy E can reach, that is <l?(r„,f) = E. 
The second equality in Equation (A2) follows from the relation 
between the delta and the Heaviside functions, that is 5(x-a) = 
^&(x-a), which also helps to simplify the calculations in Sec¬ 
tion 2.2. 

Let {X{E)) be the isotropic average of an arbitrary quantity 
Xfr, v) over the phase-space volume of a microcanonical ensemble 
of particles with energy E. From Equation (A2) this can be written 


{XiE)) = ~ f e(£-H)2f(r,v)d'H 


(A3) 


- 5(E-H)X(r,\)(fn= / /£(r,v,f)A(r,v)d'’H. 


where is the so-called the microcanonical distribution function 


/£(r,v,r)d®H = - X S{E-H)dL‘Q., (A4) 

OJ 


which defines the probability of finding a particle with an energy E 
in the phase-space volume d®H centred at the coordinates (r,v) at 
the time t (e.g. Landau & Lifshitz 1980). 


(X) = J /(r,v,f)A(r,v)d'f? (A5) 

= J J e(E-H)X(r,y)<fQ(E) 

= JJ5(E-H)X(r,\)(fn 
= J dEN(E,t) J fE(r,y)X(r,y)d'^Q 

= J N(E,t){X{E)}dE=X, 


highlighting the equivalence between phase-space average, (X), 
and average over microcanonical ensembles, X. 

For spherical & anisotropic systems the distribution function 
can be written as / = f(E,L, t), where L is the modulus of the angu¬ 
lar momentum vector. To calculate the volume accessible to stars 
with energy E and angular momentum L it is useful to define a 
coordinate system such that the tangential and radial velocity com¬ 
ponents can be written as v, = vsinS and Vy = vcosO. In these coor¬ 
dinates the angular momentum becomes L = rv sin 6 and dL = rv^d^, 
so that sin 6d9 = LdL/{rvVr). The sign of the radial velocity compo¬ 
nent changes as the particle moves toward and away from the centre 
of the potential, becoming zero at the peri and apocentre radii, rp 
and ra, respectively. With this information in mind the density of 
states of an anisotropic system can be written as 


uj(E,L,t) = / 5(E-H)S(L-rvsme)d^n 


= (47r)M r^drx v^dvsinMe 


2 f dr 2 
(47r)^L / — = SnLP(E,L,t), 

.L. Vy 


(A6) 


where P(E,L,t) = 2 dr/vy is the period of an orbit with energy 
E and angular momentum L at the time t. 

For anisotropic systems Equation (A5) becomes 


(X) = / f(r,y,t)X(r,y)d^n 


(A7) 


dEdL- / S(E-H)5(L-rvsme)X(r.y)d^n 


N{E,L, t){X(E,L))dEdL = X, 


where {X(E,L)) denotes the average of 2f for all stars of a particular 
E andL. 


APPENDIX B: MARKOV CHAINS 

Equation (30) shows that pfE, t\Eo,to) plays the role of a transition 
probability between the states f = fo and f = foH-r', where t' is a time 
interval of arbitrary length. As such, pc only depends on the initial 
and final states of the system and not on the steps taken to get there, 
a remarkable property which is worth discussing in some depth. 

Consider an evolutionary path that goes from to to to + t' 
through an intermediate step fo < ti < fo + T^ How would the 
results obtained in §3.2 change? As a first step we notice that 
since N(I, {a}^) is invariant each individual state must obey Equa- 
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tion (26), hence 


J pU\Eo,to)N(Eo,{a},,to)dEo = j p{l\Euh)N{Eu{a},,U)m (Bl) 

= J p(I\Ejo+T')N(E,{a},,to + T')dE. 
Following the same steps that led to Equation (30) we find that 

N(E,{a},,t) = J dEiN(Eu{a}.,ti)p,(E,t\Ei,h) (B2) 

= [dEoN(EQ,{a}„,to) [dEipc(E,t\Ei,ti)pc(Ei,ti\Eo,tQ). 


[f(b)+f(a)]/2, such that Equation (B7) reduces to (33). Therefore, 
we find that in potentials that evolve in a linear regime, where the 
time scale (/i//io)“* is very long compared with the orbital period of 
the particles, the transition (or “jump”) between to and to + r' given 
by Equation (30) is transitive, i.e. it neither depends on the manner 
in which the intermediate evolution proceeds nor on the length of 
r'. It also becomes clear that in order to preserve the transitivity of 
Equation (30) in potentials that vary rapidly with time the length of 
the jumps must be decreased such that the condition r' (p/po) ' 
is everywhere and at all times satisfied. 


Hence, by analogy with Equation (30) one can define a transition 
probability Pc, 2 (£,t|£'o,to) = / dEip^E,t\Ei,ti)pc(Ei,ti\Eo,to). 

The above result can be easily generalized to a serial sequence 
of K time-steps to <ti <,...,< t^-i < , with = to + In doing 

this we find that the transition probability between to and to + r' can 
be expressed as a Markov chain (e.g. Krapivsky et al. 2010) 


Pc,K(E,t\Eo,to) = J dEipc(Ei,ti\Eo,to) J dE2Pc(E2,t2\Ei,t[}... (B3) 

^ f dEfi—\Pc(^E ,, t)^-[)pc{E\\ Ei ^—2 , ti^^2^- 


We shall say that a Markov chain is transitive when the transition 
probability between the states to and to + t' is independent of the 
number of intermediate steps, that is if pc = Pc,i = ... = Pc,K for any 
value of K. 

As discussed in §3.2, transition probabilities can be ap¬ 
proached by Gaussian functions in potentials that evolve in a linear 
regime. It is straightforward to show that the mean (AC^) and the 
variance (D{) of the Gaussian resulting from the integral over j dEj 
in (B3) can be written as 

AC^ — ACfEy+i,tj+i,Ej.tj)-!-AC(ifj,,ty-i)^ (B4) 
and 

— bc{Ej+\ ,tj+i^Ejpj) + bc(Ej,tj^Ej-ipj-i) ^ 

Applying Equations (B4) and (B5) recursively from / = k - 1 to j = 
1, denoting R = R(to + T') and Ro = Rito), and using Equations (32) 
and (33) yields 

AC(E,t,Eo,to) = C(E,t)-C{Eo,to){^'^ , (B6) 


and 


D,(E,t,Eo,to) = D(E,t) + D(Eo,to)(^^'^ +2^D(Ej,tj)(^^) 


2 

t' 


D(E,x) 


R{x) 


dx, 


where the last step follows from changing ^ j dt and ap¬ 

plying the extended trapezoidal rule for k ^ 1. Equations (B6) 
and (B7) respectively show that the drift coefficient AC retains no 
memory of the evolutionary path of the system, while Dc is re¬ 
lated to a time-average of the diffusion coefficient over the interval 

topO+T . 

Intuitively, we can see that the probability pc becomes tran¬ 
sitive if the initial and final states do not deviate strongly from 
the time-averaged state of the system. This is clearly the case 
for time intervals t' (/i/po) *, for which (b — ay' f{t)dt ~ 



